Biostar：课程15、16

翻译：生物女博士

注：本系列课程翻译已获得原作者授权。

#为作者的注释

#*为译者的注释

Lecture 15 - BAM文件和samtools

# 安装一个短read模拟器。Heng Li写的wgsim。

cd ~/src

git clone https://github.com/lh3/wgsim.git

cd wgsim

gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm

ln -s ~/src/wgsim/wgsim ~/bin/wgsim

# 检查一下是否可以用。

wgsim

# 现在，下载和安装samtools。

# 把samtools链接到~/bin目录。

cd ~/src

curl -OL http://sourceforge.net/projects/samtools/files/samtools/1.1/samtools-1.1.tar.bz2

tar jxvf samtools-1.1.tar.bz2

cd samtools-1.1

make

ln -s ~/src/samtools-1.1/samtools ~/bin/

# 检查一下是否可以用。

samtools

# 查看samtools的手册。

man ~/src/samtools-1.1/samtools.1

# 我们将会利用一个突变的参考基因组来生成一些短reads，并map（回贴？）回去。

mkdir ~/edu/lec15

cd ~/edu/lec15

# 生成一些数据，并把输出保存到一个文件。

wgsim -N 10000 ~/refs/852/ebola-1999.fa read1.fq read2.fq > mutations.txt

# 运行比对。

bwa mem ~/refs/852/ebola-1999.fa read1.fq read2.fq > results.sam

# 开始转换成BAM文件。

#*运行的时候请注意，这里用的samtools版本比较旧，新版本有的命令已经改变了。

samtools view -Sb results.sam > temp.bam

# 对比对结果进行排序。

samtools sort -f temp.bam results.bam

# 索引比对结果。

samtools index results.bam

#*做完了这一步，如果想方便查看一下比对结果，可以下载一个IGV软件，并将参考序列、bam文件和bam文件的索引后生成的一个文件导入IGV即可。

# 通常，我们会把所有的命令放到一个脚本中并运行。

# 查看文末的脚本。

#*原文并未将这个脚本放上来。但是后边的课程中会有更加完善的脚本。

#*本行命令暂时不运行：bash align.sh read1.fq read2.fq results.bam

# 使用samtools过滤。

# -f 匹配标识（flag）（保留匹配的标识对应的reads）

# -F 过滤标识（flag）（去除匹配标识对应的reads，保留剩余的）

#*在上两个lecture中有说过sam格式，sam每一列有不同含义，其中一列（第二列）叫标识（flag），不同数字有不同含义

#*下边内容引用自参考资料1

FLAG, 概括出一个合适的标记，各个数字分别代表

1     序列是一对序列中的一个

2     比对结果是一个pair-end比对的末端

4     没有找到位点

8     这个序列是pair中的一个但是没有找到位点

16   在这个比对上的位点，序列与参考序列反向互补

32   这个序列在pair-end中的的mate序列与参考序列反响互补

64   序列是 mate 1

128 序列是 mate 2

假如说标记为以上列举出的数目，就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字，比如说83=（64+16+2+1），就是这几种情况值和。

# 比对到反向互补链的reads。

samtools view -f 16 results.bam

# 比对到正向链的reads。

samtools view -F 16 results.bam

# -c 可以用来计数

samtools view -F 16 results.bam | wc -l

samtools view -c -F 16 results.bam

# 用质量来过滤。BWA 的mapping质量为0时，表示reads是map到多个位置。而 q>=1表示该reads是单一mapping。

samtools view -c -q 1 results.bam

# 统计高质量的比对数目。

samtools view -c -q 40 results.bam

# 参考基因组的名字为 gi|10313991|ref|NC_002549.1|

# 总是手动输入会有些乏味，我们可以创建一个简称。

CHR='gi|10313991|ref|NC_002549.1|'

# 截取部分数据。

samtools view results.bam $CHR:1-100

# 查看文件的特定某列

samtools view results.bam $CHR:1-100 | cut -f 4 |  more

Lecture 16 - Converting data with ReadSeq

# 安装文件格式转换器。

# 为之创建一个目录。

mkdir -p x

cd ~/src/readseq

curl -OL http://iubio.bio.indiana.edu/soft/molbio/readseq/java/readseq.jar

# 跟调用trimmomatic很相似

# 将调用trimmomatic脚本修改一下。

#*在lecture8时，我们创建了一个文件叫trimmomatic来简化trimmomatic的命令执行。详情参见Biostar：课程7、8

cp ~/bin/trimmomatic ~/bin/readseq

# 用下边程序来替换：

# java -jar ~/src/readseq/readseq.jar $@

# 你也可以加上：

# java -cp ~/src/readseq/readseq.jar app $@

# 用图形用户界面来运行。

# 获取完整的GenBank记录的1999年的埃博拉基因组。

# http://www.ncbi.nlm.nih.gov/genome/4887

# 获取完整genbank文件。

mkdir -p ~/edu/lec16

cd -p ~/edu/lec16

# 获取完整数据。

#*efetch貌似有问题，可以用网页直接下载。详情参见Biostar：课程3、4

efetch -db nucleotide -id NC_002549.1 -format gb > NC.gb

# 提取数据未GFF（Generic Feature Format）格式。

# 自动获取输入文件的格式。

readseq -format=GFF NC.gb

# 你也可以设置输出文件的名字。

readseq -format=GFF -o NC-all.gff NC.gb

# 提取为fasta文件。

readseq -format=FASTA -o NC.fa NC.gb

# 从GFF文件中提取“gene”的行

#* \t表示的是tab分隔符

cat NC-all.gff | egrep '\tgene\t'

# 同时，前两行也是我们要的。

echo '##gff-version 2' > NC-genes.gff

cat NC-all.gff | egrep '\tgene\t' >> NC-genes.gff

# 染色体坐标不匹配！我们需要去把它修复一下。

# 我们需要重做比对，或者交换现有比对或特征的坐标。

# 重新索引和比对数据。

cp NC.fa ~/refs/852/

bwa index ~/refs/852/NC.fa

cp ../lec15/*.fq .

# 编辑align.sh使其能用于新文件，并重新运行。

bash align.sh read1.fq read2.fq results.bam

# 把数据导入IGV并查看。

# 查看100-500bp区域的深度。

# 输出为序列号、基因组索引和覆盖度。

samtools depth -r NC_002549:100-500 results.bam

参考资料1：http://www.cnblogs.com

Biostar非常适合有生物学基础且想自学生信的这部分人入门。

Biostar课程文章系列：

欢迎关注我们